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The collision density of photons in an infinite Compton-scattering medium was calcu- 
lated by random sampling. The intensity, spectral composition, and angular distribution 
were then determined for radiation reflected and transmitted by plane-parallel barriers by 
carrying out appropriate integrations over the collision density, taking into account the 
absorbing properties of the medium and the effect of boundaries. All numerical work was 
done on the National Bureau of Standards automatic computer (SEAC). Sample results 
are presented for 0.66-Mev radiation incident on water barriers of various thicknesses. 



1. Introduction 

The Problem. The results of an investigation of 
the following boundary problem in gamma-ray 
diffusion are presented: A beam of monoenergetic 
gamma rays is incident at a given angle on a plane- 
parallel barrier that has a finite thickness in one 
dimension but is infinite in the other two dimensions. 
What are the intensities and spectral compositions 
of the reflected and transmitted beam? 

Method of Solution. An analytical attack on this 
problem by the solution of the relevant transport 
equation leads to great mathematical difficulty. 
So far such an approach has been successfully carried 
out for homogenous infinitely extended media [l], 3 
but boundary problems remain unsolved. In order 
to bypass these difficulties, explore the problem, and 
obtain at least an approximate solution, this investi- 
gation employed the Monte Carlo (random sampling) 
method. All numerical work was carried out on the 
automatic computer of the National Bureau of Stand- 
ards (SEAC). 

The calculation was divided, in regard to both 
method and execution, into two rather distinct parts: 
(1) a stochastic part, in which photon random walks 
are generated by random sampling in an infinite 
nonabsorbing scatterer, whereby a collision density 
is obtained, and (2) a calculation of the reflected 
and transmitted radiation obtained by summations 
over collision densities, taking into account boundary 
conditions and absorption characteristics of the 
barrier. There were a number of reasons for this 
division. For one thing, it reduced the required 
computer memory capacity. More importantly, 
the random walks have a "universal" character, 
i. e., they can be used for the calculation of different 
boundary problems involving different geometries 
as well as different scattering and absorbing media. 
Not only will the repeated use of the same set of 

i This work was supported by the Office of Naval Research and the Atomic 
Energy Commission Reactor Division. 

2 Some of the results of this paper were presented at a symposium on Monte 
Carlo methods at the University of Florida in March 1954. 

^Figures in brackets indicate the literature references at the end of this paper. 



random walks often lead to computing economy, 
but it may also increase the accuracy of a calculation. 
Suppose, for example, that we wish to determine the 
difference in transmission for two beams incident on 
a barrier at different angles but otherwise identical. 
It will then be to our advantage to base the compari- 
son on the same set of random walks because 
irrelevant differences resulting from statistical fluc- 
tuations will thereby be largely eliminated. 



2. Gamma-Ray Random Walks in an 
finite Compton Scatterer 

2.1. Definitions 



In- 



The state S of a photon can be specified by a set 
of six quantities: 

S=(E,e,<p,x,y } z), 

where E is the photon energy, and (p are angular 
coordinates describing the direction of motion (in a 
spherical coordinate system with the 2-axis as polar 
axis), and x, y, and z are Cartesian coordinates of 
position. Let the state of the photon immediately 
after the nXh scattering event occurring in a given 
random walk be denoted by S n (n=l,2,- • •)> and 
let So denote the state in which it was introduced 
into the scattering medium. A random walk is 
then described by a sequence 



Of), Oi, O2, 



& 



each term (except So) depending stochastically on 
its immediate predecessor only (Markov process). 
The length, L, of such a sequence would be infinite, 
but in the present work the random walks are ter- 
minated when the energy, E } drops below 30 kev. 
This arbitrary cutoff 4 results in an average value 
<L>=18 for an initial energy #0=660 kev. A 



4 The cutoff is justified ou physical grounds because radiation at energies below 
30 kev is always so heavily absorbed that it makes only a negligible contribution 
to the emergent radiation flux. 
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set of random walks constitutes a Monte Carlo 
estimate of the collision density resulting from the 
diffusion of radiation in an infinite homogenous 
Compton scatterer. 

It will be observed that for the problem of trans- 
mission and reflection by plane-parallel barriers, 
only the variables E, 6 and one space variable, say z, 
are required. Nevertheless, it is worthwhile to 
calculate the other three variables also. The ran- 
dom walks are thus made applicable to other 
boundary problems with different geometrical condi- 
tions. Moreover, we can use the same set of 
random walks for barrier problems with different 
angles of incidence. We always set 6 =0 and rotate 
the boundaries as required b}^ using, instead of cos 
6 and z, 

cos 0'=cos 6 cos a+sin 6 sin a cos <p 

z' = x sin a-\-z cos a, 

where a = 9' is the angle of incidence. 

2.2. Random Sampling 

Next the sampling scheme is described, by which 
successive states S n are selected, the initial state 
being specified as $0= (E Q , 0,0, 0,0,0). Calculations of 
a set of random walks with this initial condition will 
serve for the solution of problems involving mono- 
energetic beams of energy E and arbitrary direction 
of incidence a. Because of the linearity of the 
gamma-ray diffusion equation, it is possible to 
obtain solutions for incident beams of specified 
energetic and angular composition by superposition 
of the results obtained with different E '& and as. 

Prior to the discussion of the detailed steps in the 
sampling process, some comment is in order concern- 
ing the required random numbers. With a high- 
speed computer the use of tables of random numbers 
is clearly impractical. Instead, so-called pseudo- 
random numbers r m (Q<r m < 1) were used, and which 
were generated as required in the course of the 
computation by a method developed by O. Taussky- 
Todd [2]. They are defined by the relations 



r m — " *+ m 

i? m+ i=5 17 i? m modulo2 42 , B =l. 



(2) 



It may be shown that the period of r m is 2 40 , i. e., 
that a sequence of 2 40 different numbers will be 
obtained before repetition occurs. Extensive test- 
ing carried out at the National Bureau of Standards 
has shown that these pseudorandom numbers 
satisfy the various accepted statistical criteria of 
randomness. It is an advantageous feature of this 
method that identical random walks may be re- 
created repeatedly for checking purposes, provided 
the initial random number used for the walk is 
recorded. 



The various steps 
given S n , follow. 



in the calculation of S, 



a. Energy Change 

The energy E rc+l after a scattering is determined by 



f £ n k(E n ,E)dEr n k(E m E)dE\=r, 



(3) 



where r is a random number, and k(E n ,E) is the 
Klein-Nishina differential coefficient (per unit path 
length) for Compton scattering with energy change 
from E n to E. 

It was found convenient to carry out the random 
sampling for the wavelength \=mc 2 /E, where mc 2 is 
the electron rest mass, rather than for the energy 
directly. If we make this change of variable, sub- 
stitute the explicit formula for the Klein-Nishina 
differential coefficient, and carry out the indicated 
integrations, eq (3) is transformed into 



[l-2X n (X B +l]log^- 



■X.(X.+2) 



X«+i A„, 



Xn(A n +l~-A M ) + 



2X 2 



H [l-2X B (X»+l)]log^- 



'(2+A„) 2 



(4) 



n+1) 



To use this equation to find X n+1 , given X^ and r, 
would be intolerably cumbersome, even when an 
automatic computer is used. Latter and Kahn [3] 
evaluated eq (4) numerically and presented a table 
of X n+ i=X w+ i(X w ,r) for a grid of values in the range 
0<r<l, 0.05<X„<10. The tabulation has been 
extended to X„ = 16 and used to obtain a polynomial 
representation 

\.+i=SS4WV. (5) 

i=0 j=0 

This required extensive bivariate interpolation, 
which was carried out by a convenient matrix tech- 
nique. This technique is described in the appendix 
because it does not appear in the standard literature 
on numerical analysis. 

The following coefficient-matrices, A ih were found 
to give representations of X ?+1 accurate to 1 percent: 



0.05<X n <0.5 


i 






J 









1 


2 


3 


4 


5 



1 
2 
3 




1.0 




0. 1180 
. 07930 
2.124 
-3. 418 


-0. 9844 
20.57 
-63.66 
70.31 


4.683 
-69. 56 
209.7 
-227. 8 


-6.896 
117.2 
-342.3 
359.2 


5.079 

-68. 29 

194.1 

-198. 3 


0.5<X n <10.0 


i 


1 
2 
3 


3 






1.0 




1 


2 


3 


4 


5 


0. 04100, 
.4779 

-.05409 
. 002153 


2.851 
-1.502 
0. 07590 
. 0008490 


-8. 346 

4.863 

-0. 1713 

-. 0077 


13.37 
-5. 168 
-0. 02945 
. 02359 


-5.916 
1.330 
0. 1789 
-. 01891 



344 



10<X»<16 


i 


./' 





1 


2 


3 



1 
2 
3 




1.0 




27. 4295 
-6. 173625 
0. 4584375 
-.01106 


-96. 637". 
23. 35 
-1. 763 
0. 04342 


71. 2080 
-17.17 
1.304 
0. 03233 



b. Change of Direction 

Deflection and wavelength shift are related by the 
Compton equation 



cos u n =l — \ n+1 + \ n , 



(6) 



where co n is the angle between the directions of motion 
immediately before and after the n+lst scattering. 
The direction of the scattered photon depends on 
o) n and on an azimuth angle x n , which is a random 
quantity distributed uniformly between and 2tt. 
We need the sine and cosine of x n rather than x n 
itself. Hence we take advantage of the following 
convenient computational scheme suggested by Von 
Neumann [4]. Choose random numbers a and b 
satisfying the condition a 2 + 6 2 <l, and let c be a 
random number that is equal to ±1 with probability 
1/2. Then 



2abc 

COS X n =-2 



a 2 -b 2 



,„ and sin x„=- 9 
a-^rb 2 a 2 +b 2 



(7) 



From the deflections co„ and x n , the new angular 
coordinates are then determined by the following 
trigonometric relationships: 



cos fl^+^cos B n cos w w +sin 6 n sin oo n cos x n 
sin x n sin co „ 



sin(<£ n+1 — <f> n )= 
cos(0„ + i — 4>„y 



sin e n+l 
cos co n — cos d v cos 6 



w+l 



r (8) 



sin B n sin n+1 

c. Displacement 

The position of the (n+ l)st scattering is defined by 



(9) 



%n-\-l %n 


sin u n <^ v n j Qtr r 

v* s (E n ) 


Vn+1 — Vn 


sin 6 n sin cp n -, 


%n+l — %n 


- C0S *»Wr 
n,(E n ) 10gr 


srhere r is a random 


number, and 



!*„(E n )=r n k(E n ,E)dE. 



3. The Transmission-Reflection Boundary 
Problem 

3.1. Summation Over the Collision Density 5 

The information obtained by sampling J photon 
random walks is contained in the set of state vectors 
S n (j)(n=0,l,- • - Lj; j =1,2,- • • J), which describe the 
state of the 7th photon immediately after its nth 
scattering. We consider these vectors to form a 
set of representative points (a sample collision 
density) in a six-dimensional collision-density space. 
The solution of boundary problems can be obtained 
by performing appropriate sums over the collision 
density. 

In the reflection-transmission problem, one con- 
siders a plane-parallel barrier between the planes 
z=0 and z=t. A beam of photons is assumed to be 
incident on the face 2=0. 

Reflected and transmitted radiation will be classi- 
fied according to energy and direction. 6 Let R ik 
denote the average number of photons per incident 
photon that are reflected from the face 2=0 with 
energies in the ^th interval and directions in the Hh 
interval. Similarly, let T ik denote the average num- 
ber of photons transmitted through the face z=t. 
E ik and T ik are obtained by summations over the 
sample collision density as follows: 



■hik — 






The factors in this formula have the following 
meaning: 

(a) The factor B n (j) takes into account the bound- 
ary conditions. It is equal to unity when the in- 
equalities 0<z m (j)<t are satisfied for m = 0,], • • • n, 
and zero otherwise. In other words, a contribution 
to the reflected radiation is obtained only if the 
photon in the state under consideration has not 
previously gone out of the barrier. 

(b) The factor 



n-l 

Pn(j)= H exp 

TO=0 



{- 



sec 6 m ( j) 

»A[E m (j)} 



[2, 



,+i a)- z m (:/)] y 



en) 



is the probability that the photon has not been 
absorbed in the jth. history prior to reaching state 
S n (j)- (ma is the absorption coefficient.) 

(c) The quantity Q n (j,z') denotes the probability 
that a photon starting from state S n (j) will cross the 
plane z=z' without further scattering or absorption. 
It is defined as follows: 

Qn(j,z') vanishes except when sec 6 n (j)[z' — z n (j)] 
>0, in which case 



Q„0>') exp 



{• 



s ec 6„(j) 
'»[E n (j)] 



z —z. 



.(/)]} 



(12) 



5 The method described in this section is closely related to suggestions by 
H. Kahn, [3]. 

6 32 energy intervals and 10 angular intervals have been used. 
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(/*= Ma+ Ms is the total linear attenuation coefficient). 

(d) /i fc (i) is the energy-angle classification factor. 
It is equal to unity when E n (j) is in the ith. energy 
interval and B n (j) in the kth. angular interval, and 
zero otherwise. 

Transmitted radiation T ik is calculated by a for- 
mula similar to eq (10) but with Q n (j,0) replaced by 

Qn(j,t). 

We shall use the collision density also for the eval- 
uation of the flux in an infinite medium in which a 
plane source is embedded in the plane 2=0. By 
comparing the results for this problem with those 
for the barrier, the nature of the boundary effect can 
be established. Moreover, the Monte Carlo results 
for the infinite medium can be checked against those 
obtained by systematic moment-method calculations. 
Such a comparison is of use in establishing the valid- 
ity and accuracy of the random-sampling approach. 
Let B ik (z f ) denote the average number of photons 
crossing the plane z—z' from both directions (our 
main interest will be in the planes that form the 
barrier boundaries: z' = and z' = t). The necessary 
summation over the collision density is similar to 
formula (10), except that the boundary factor 
B n (j) is omitted. 

When the sampled collision density is used for the 
simultaneous determination of R ik , T ik , B ik (0), and 
B ik (f), all states S n (j) make a contribution to the 
final score except those for which z n (j)<^0 and 
cos0 w (j)<O, or z n (j)y>t and cos w (j)>O. This is in 
contrast to a direct stochastic analogue method, 
which would allow only those sections of the random 
walks containing an actual crossing of the boundaries 
to contribute to the score. The more elaborate 
scoring procedure gives proper credit for "near 
misses/' i. e., collisions taking place close to a 
boundary, for which the factor Q n (j,z') is very close 
to unity. 

3.2. Fluxes and Buildup Factors 

The SEAC calculation is designed to yield the 
quantities R ik , T ik , and B ik , which represent the 
average number of photons crossing the planes 2=0 
and z=t per unit area (the number flux). Radi- 
ation detectors often do not measure the number 
flux directly, so that it is desirable to calculate 
various related types of flux. These can all be 
obtained from the number flux through multi- 
plication by suitable weight factors, depending on 
the energy and direction of the photons. For 
example, multiplication of the number flux by the 
average energy, (E t ), yields the energy flux, division 
by the average value, absolute value |{cos d k }\, 
yields fluxes through a unit area whose normal is in 
the direction of motion of the radiation. This is 
the flux that a nondirectional detector would 
measure. Some inaccuracy is introduced into these 
flux-conversions because the averages (E t ) and 
(cos k ) are to some extent arbitrary. Arithmetic 
means have been used. 



According to the usage common in shielding 
calculations the results for the total flux (integrated 
over-all directions and spectral energies) are pre- 
sented in the form of buildup factors, defined as 
the ratio of the total to the unscattered flux [5]. 
For the various types of flux described above, the 
buildup factors for an infinite homogenous medium 
have the following form: 

Number buildup factors 

B N =^2B ik /exp (— mo^/cos a) 

ik 

-nt = y^ B ik / exp (— vqz/cos a) , . 

N ir\(cosd k )\l cos a 

Energy buildup factors 
B E =QO l {E i )B ik /E exp (-^z/cos a) 

ik 



ik I (cos k ) | / cos a 



exp (— /x z/cos a), (14) 



where /x is the total attenuation coefficient of the 
source radiation, and z is the distance between the 
source plane and the plane of observation. In 
order to obtain the corresponding buildup factors 
T N , T f N , T E , and T E , one replaces B ik in eq (13) 
and (14) by T ik , and z by the barrier thickness, t. 
The reflection buildup factors R N , R f N , R Ej and R' E 
are obtained by using R ik and setting 2=0. 

4. Results 7 

4.1. Energy Spectra 

In figures 1, a, b, c, and d, energy spectra are 
shown in histogram form for the number flux of 
photons reflected and transmitted by barriers, 
reflected by a semi-infinite medium, and for the 
flux in an infinite medium through planes that 
correspond to the barrier boundary planes. The 
histograms pertain to scattered radiation only. 8 
Spectra are shown for normal incidences (a=0°) 
and for oblique incidence (a=60°). The oblique bar- 
rier thickness is the same for both cases : (/zrf/cos a) = 2. 
The shaded areas between the "finite slab" and the 
"infinite medium" histograms are a measure of the 
boundary effect. It can be seen that the boundaries 
greatly reduce the amount of soft radiation below 
200 kev. The reason for this is that in an infinite 
medium soft radiation is mainly due to photons that 
have overshot the boundary plane and have been 
turned around in a collision in the close vicinity of 
the boundary, making their second passage through 
the boundary with resultant low energy. 

7 Based on the analysis of 400 random walks for each problem. 

8 To obtain the scattered flux only, one must begin the sum in eq (10) with 
n=l. 



346 



.18 



.16 



o 

H- 
O 

I .141- 

| .12- 

o 

z 

a. 

<j) 

§ .08- 

»- 

o 

°~ .06- 

Q 
Ul 

a .04- 

_j 

ui 

01 .02- 



1 


1 1 1 1 1 1 1 
E =660 kev 


I 


" 


a - 0° 
cos a 






WATER 




; 


i 


^ 


INFINITE MEDIUM 


- 


- 




L 


SEMI- INFINITE MEDIUM 


- 


- 










yr- FINITE SLAB 


- 


- 












— 1 ° 


- 




VA 








! 


i i 


\— i i i 


] 



i r 



INFINITE MEDIUM 




"i 1 1 r 

E = 660 kev 
a ■ 60° 

cos a 
WATER 



SEMI- INFINITE MEDIUM 



.FINITE SLAB 



^ 



J L 




JZL 



z .16- 



.08- 



o .06 



.04- 



.02 



1 1 1 1 1 i r 




"I T 



1 1 1 

E - 660 kev 
a ■ 60° 

t£±_ , 
cos a 

WATER 



420 480 540 600 




. kev 




420 480 



600 



Figure 1. Energy spectra of photons reflected and transmitted by water barriers and of photons in an iniinite homogenous medium 

crossing planes corresponding to the barrier boundaries. 

The spectra pertain to a monodirectional beam of 0.66-Mev photons incident on the barrier (or released within a semi-infinite medium) with obliquity a. 

a, Reflection; obliquity a=0° (normal incidence); b, Reflection; obliquity a=60°; c, Transmission; obliquity «=0° (normal incidence); d, Transmission; obliquity 

a = 60°. 



4.2. Angular Distributions 

Figure 2, a, shows the angular distribution of the 
reflected number flux from a barrier with oblique 
thickness /zo£=l, and from a semi-infinite medium 
for obliquity angles of incidence a=0° and a = 60°. 
Figure 2, b, shows the angular spectra of transmitted 
scattered photons for normal incidence and barrier 
thickness fJtot=l and /x £=4. 



4.3. Buildup Factors 

The various buildup factors for transmission, 
reflection, and an infinite medium (as defined in 
section 3.1 and 3.2) are listed in table 1 for the 
number flux, and in table 2 for the energy flux for 
0.66-Mev radiation in water. Estimated standard 
deviations are indicated, which were obtained by 
dividing the 400 photon histories into groups of 50 
each and computing the dispersion of the buildup 
factors obtained for the various groups. 
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Figure 2. Angular distributions of photons reflected and transmitted by water barriers. 

The distributions pertain to a raonodirectional beam of 0.66-Mev photons incident on the barrier with obliquity a. 
a, Reflection; obliquities a=0° and 60°; b, Transmission; obliquity a=0° (normal incidence). 



Table 1. Number buildup factors for 0.66-Mev radiation in 

water 





not 














a 




Rn 


K'n 


T N 


T ' K 


B N 


*'* 






COS a 













































1.61 
±0.09 


2.30 
±0.16 





1 


1.18 


1.36 


2.09 


2.73 


3.69 


6.24 






±0.03 


±0.05 


±0.16 


±0.18 


±0.26 


±0.35 





2 


1.29 


1.55 


3.94 


6.02 


7.11 


12.3 






±0.03 


±0.06 


±0.45 


±0.55 


±0.64 


±1.3 





4 


1.31 


1.59 


8.21 


13.5 


13.3 


24.3 






±0.05 


±0.07 


±1.75 


±2.2 


±2.8 


±4.2 





CO 


1.31 
±0.05 


1.61 
±0.07 










60 













1.85 
±0.11 


1.97 
±0. 13 


60 


1 


1.27 


1.30 


2.00 


2.17 


4.13 


4.41 






±0.05 


±0.06 


±0.13 


±0.13 


±0.25 


±0.28 


60 


2 


1.41 


1.48 


4.13 


3.90 


8.89 


8.21 






±0.06 


±0.07 


±0. 46 


±0.50 


±0. 95 


±0.98 


60 


4 


1.47 


1.53 


15.4 


14.2 


29.3 


30.2 






±0.06 


±0.07 


±2.0 


±2.5 


±4.9 


±5.3 


60 


CO 


1.48 
±0.06 


1.53 

±0.07 











Table 2. 



Energy buildup factors for 0.66-Mev radiation in 
water 





IJ.4 


Re 


R B 


T E 


n 


Be 


K 








COS a. 














o 















1.10 
±0.02 


1.23 
±0.04 





1 


1.05 


1.10 


1.62 


1.89 


1.87 


2.48 






±0.01 


±0.02 


±0.06 


±0.09 


±0.07 


±0.10 





2 


1.07 


1.14 


2.50 


3.21 


2.97 


4.07 






±0.02 


±0.03 


±0.20 


±0.30 


±0.22 


±0.31 





4 


1.07 


1.15 


3.62 


5.64 


4.33 


6.98 






±0.02 


±0.03 


±0.39 


±0.65 


±0.39 


±0.72 





CO 


1.07 
±0.02 


1.15 
±0.03 










60 













1.19 
±0.03 


1.23 
±0.04 


60 


1 


1.10 


1.11 


1.57 


1.65 


1.98 


2.15 






±0.03 


±0.05 


±0.07 


±0.09 


±0.09 


±0.10 


60 


2 


1.13 


1.17 


2.43 


2.37 


3.24 


3.48 






±0.03 


±0.05 


±0.22 


±0.20 


±0.25 


±0.29 


60 


4 


1.14 


1.18 


5.92 


5.40 


7.61 


6.21 






±0.03 


±0-05 


±0.68 


±0.73 


±0.91 


±1.05 


60 


CO 


1.14 
±0.03 


1.18 
±0.05 
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5. Discussion 

The results presented have a preliminary charac- 
ter. Greater sample sizes are desirable, as well as 
the extension of the calculation to other energies 
and media. Moreover, the scoring method used is 
not expected to be very successful for penetrations 
of depth mo#^>4 because excessively large sample 
sizes would be required. Prior to extending the 
calculations by the collision-density method, other 
Monte Carlo approaches have been explored, in- 
cluding a semianalytic method designed for dealing 
with deep as well as shallow penetrations [6]. Table 
3 gives a comparison of the buildup factors obtained 
in this investigation (a) with those found by the semi- 
analytical Monte Carlo method, (b) with the Monte 
Carlo results of a group at the Naval Research 
Laboratory [7], and (c), for an infinite medium, with 
calculations carried out according to the moment 
method of Spencer and Fano [1, 8]. The agreement 
is generally good, indicating that the shallow 
penetration of gamma radiation, in the presence of 
boundaries, can be calculated accurately by the 

Table 3. Comparison of buildup factors obtained in various 
calculations for 0.66-Mev radiation in water 

Column 1, This paper; 2, semianalytic Monte Carlo method [6]; 3, Monte 
Carlo, NRL [7]; 4, Spencer-Fano moment method [8]. 
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Monte Carlo met luxl. The results for R N and B E 
at 0.66 Mcv are consistent with the results of a 
Monte Carlo calculation at 1 Mev by Hayward and 
Hubbell [9]. 
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7. Appendix. A Method of Bivariate Inter- 
polation 

Suppose we know a bivariate function f(j>,q) a1 
(M+1)(N+1) gridpoints (j) St q t ) (s=0,l, v • , M; 
£=0,1, . . . , N). The problem is to use this knowl- 
edge for determining the expansion coefficients in a 
polynomial expansion 

/(P,5)=SSV2'. (17) 

i = j = 

Let F=f(p 8} q t ) be the known matrix of functional 
values, and A=(Aij) the desired coefficient matrix. 
From the variables p and q we form the matrices 
P=(psY and Q=(qtY. In matrix form, eq (17) 
becomes 

F=PAQ', (18) 



where Q' is the transpose of Q. Hence 
A=P- l F(Q- 1 )'. 



(19) 



P and Q are so-called alternant matrices whoso in- 
verse is well known [10]. The prescription for ob- 
taining' the inverse is simplest to describe by means of 
an example. When P has rank 3 



r l Po p^ 

1 Pi P\ 

<X P2 pV 



P1P2 



Pop?, 



P0P1 



(VO — P1)(PO — P>2) (Pl — po)(Pl — V2) (P2-PO)(P2 — PO 

~(Pl + P2) (P0 + P2) —(P0+Pl) 

(P0 — P1)(P0 — P2) (P1 — P0)(PI~P2) (P2~P0)(P2 — PI) 

1 1 1 

L(P0 — Pi)(P0 — P2) (Pl — po)(Pl — V2) (P2 — po)(P2 — Pl). 



(20) 
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The numerators of of the elements in the respective 
columns are elementary symmetric functions in the 
arguments p Q , p x , and p 2 , with one argument omitted 
each time. The arrangement of the denominators 
is obvious. 

When the grid points are equidistant, one can, 
without loss of generality, set p s =s. The inverses 
P _1 for this case are listed below for the first few 
alternant matrices : 



The sum of the elements in the first row of each 
inverse matrix is 1, and the sum of the elements in 
each other row is zero. This is also true for the case 
of nonequidistant grid points, and provides a useful 
arithmetical check. 



The author thanks Frank Stockmal for his help in 
coding the problem for SEAC, and Mary Orr and 
John Doggett for their assistance with the hand 
computations. 
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